function [Final_V,B,dnormarray,S]=DNMFCSLT(X,U,V,S,B,layers,options)

r1=options.eta1;

lambda1=options.lambda1;
lambda2=options.lambda2;
lambda3=options.lambda3;
cluster_num=options.cluster_num;
mu=1e-4;
rho=1.6;
mu_max=1e10;
epsilon=options.tolfun;
maxiter = options.maxiter;
[m_fea,n]=size(X);
m = numel(layers);
F=zeros(layers(m),cluster_num); 
tau = ones(m, 1) / m;
Q = cell(1, m + 1);
tmpDist = cell(1, m);
% S=cell(1,m);
R=cell(1,m);
G=cell(1,m);
A=cell(1,m);
D=cell(1,m);
L=cell(1,m);
beta=cell(1,m);
Z=cell(1,m);
M=cell(1,m);

for v=1:m
	R{1,v} = zeros(n,n);
    G{1,v} = zeros(n,n);
end


dnormarray = [];
dnormarray(1) = realmax;
iter = 1;
Isconverg = 0;
while(iter<=maxiter)
	fprintf('----processing iter %d--------\n', iter);
	Q{m + 1} = eye(layers(m));
	for i_layer = m:-1:2
        Q{i_layer} = U{i_layer} * Q{i_layer + 1}; 
    	end
	
	for i=1:m
		A{i}=(S{i}+S{i}')./2;
		
		D{i} = diag(sum(A{i}));
		L{i}=D{i}-A{i};
	end
	
	%%VmVmT = V{m} * V{m}'
	
	for i=1:m
		%%更新Ui
		if i==1
			Ru=X * V{m}' * Q{2}'; 
			Rd=U{1} * Q{2} * V{m}*V{m}' * Q{2}'; 
			U{1} = U{1} .* Ru ./ max(Rd, 1e-10); 
		else
			Ru=P' * X * V{m}' * Q{i+1}'; 
			Rd=P' * P * U{i} * Q{i+1} * V{m}*V{m}' * Q{i+1}';
			U{i} = U{i} .* Ru ./ max(Rd, 1e-10); 
		end
		
		if i == 1
            		P = U{i}; 
        	else 
            		P = P * U{i};
       		end
       		
       		for i_layer = m:-1:2
        		Q{i_layer} = U{i_layer} * Q{i_layer + 1}; 
    		end
       		
		
		%%更新Vi
		if(i<m)
		Vu=P' * X + lambda1*tau(i)^r1 * V{i} * A{i}; 
		Vd=P' * P * V{i} + lambda1*tau(i)^r1 * V{i} * D{i}; 
		V{i}=V{i} .* Vu ./ max(Vd,1e-10);
		else
		%%更新Vm
		Vu=P' * X + lambda1*tau(i)^r1 * V{i} * A{i}+lambda3*F*B; 
		Vd=P' * P * V{i} + lambda1*tau(i)^r1 * V{i} * D{i}+lambda3*V{i}; 
		V{i}=V{i} .* Vu ./ max(Vd,1e-10);
		end
		%%更新Si
	
		distx = [];
		distx = EuDist2(V{i}',V{i}',0);
		distx=distx-diag(diag(distx));
		beta{i}=distx;
		beta{i}=beta{i}.*(lambda1*0.5*tau(i)^r1);
		Z{i}=G{i}-R{i}./mu;
		M{i}=Z{i}-beta{i}./mu;
		S{i} = zeros(n);
		
		temp_j=1:n;
		for j=1:n
			S{i}(j,temp_j~=j) =opt_S(M{i}(j,temp_j~=j));
		end
	end
	%%更新权重
	for v = 1 : m
		tmpDist{v} = EuDist2(V{v}',V{v}',0);
	end
	tmpTau = tau;
	for v = 1 : m
		tmpTau(v) = sum(sum(tmpDist{v} .* S{v}));
	end
	tau = (tmpTau + eps).^(1 / (1 - r1));
	tau = tau / (sum(tau));
	
	% Update F
	nc = sum(B, 2); 
	F=bsxfun(@rdivide,  V{m} * B', nc'); 
	emptyClusters = (sum(B') == 0);
	F(:, emptyClusters) = 0;
	
	
	% Update B
	norms = sum(F.^2);
	distance = bsxfun(@plus, -2 * F' * V{m}, norms'); 
	[distance, labels] = min(distance);
	B = full(sparse(labels, 1 : n, 1, cluster_num, n, n));
	
	
	
	
	%%更新G
	S_tensor = cat(3, S{:,:});
    R_tensor = cat(3, R{:,:});
    temp_S = S_tensor(:);
    temp_R = R_tensor(:);
	sX = [n, n, m];
    %twist-version
    [g, objV] = Gshrink(temp_S + 1/mu*temp_R,(n*lambda2)/mu,sX,0,3)   ; %%%%%%%%
    G_tensor = reshape(g, sX);
	
	%%更新R
	temp_R = temp_R + mu*(temp_S - g);
	R_tensor = reshape(temp_R , sX);
	
	
	%%更新mu
	mu = min(rho*mu, mu_max);
	
	for v=1:m
		G{v} = G_tensor(:,:,v);
		R{v} = R_tensor(:,:,v);
	end
	
	
	
	dnormarray(iter)=calculate_loss(X,U,V{m},L,lambda1,m,tau,F,B,lambda3);
	if(dnormarray(iter)<epsilon)
		break;
	end
	iter = iter + 1;
	Final_V=V{m};
end
end

function loss=calculate_loss(X,U,Vm,L,lambda1,m,tau,F,B,lambda3)
	P=U{1};
	for i=2:m
		P=P*U{i};
	end
	loss=norm(X -  P* Vm, 'fro')^2;%+lambda1/2*loss_graph+lambda3*norm(Vm-F*B,'fro')^2;
	loss=sqrt(loss);
end
